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ABSTRACT 

The  reliability  polynomial  of  a  graph  gives  the  probability  that  a  graph  is  connected  as  a  function  of  the 
probability  that  each  edge  is  connected.  The  coefficients  of  the  reliability  polynomial  count  the  number 
of  connected  subgraphs  of  various  sizes.  Algorithms  based  on  sequential  importance  sampling  (SIS)  have 
been  proposed  to  estimate  a  graph’s  reliability  polynomial.  We  develop  a  new  bottom-up  SIS  algorithm 
for  estimating  the  reliability  polynomial  by  choosing  a  spanning  tree  and  adding  edges.  This  algorithm 
improves  on  existing  bottom-up  algorithms  in  that  it  has  lower  complexity  ~  CHE1)  as  opposed  to  0(EV3), 
and  it  uses  importance  sampling  to  reduce  variance. 

1  INTRODUCTION 

Let  G  be  a  connected  graph  with  vertex  set  V  and  edge  set  E.  We  define  r(x),  the  reliability  polynomial 
of  G,  to  be  the  probability  that  the  graph  remains  connected  when  edges  arc  removed  independently  with 
probability  x.  This  function  is  a  polynomial  which  can  be  written  as 

|£|-|V|+1 

r(x)  =  Y,  Njxfl-x)1 
i= o 

where  /V,  is  the  number  of  connected  subgraphs  of  G  with  \E\  —  i  edges. 

It  will  be  more  convenient  for  us  to  work  with  this  factored  form.  We  define  the  reliability  generator 

p(x)  =  YNiX'- 

i 

In  particular,  we  use  the  term,  “low”  order  coefficient  to  refer  to  Nj  for  small  i,  that  is,  the  number  of  ways 
of  removing  just  a  few  edges  from  G  and  still  have  the  resulting  graph  connected.  And  we  use  “high” 
order  coefficient  to  refer  to  Nj  for  large  /,  that  is  the  number  of  connected  graphs  that  can  be  formed  by 
adding  a  few  edges  to  a  spanning  tree  of  G.  The  reliability  generator  can  be  used  to  evaluate  the  reliability 
polynomial  itself  via 
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Computing  the  reliability  of  a  graph  is  computationally  hard  (Ball  1986).  A  variety  of  algorithms  have 
been  proposed  to  approximate  the  reliability  of  a  graph.  Some  of  these  algorithms  seek  to  compute  r(x o) 
or  1  —  r(xo)  for  some  fixed  probability  xo,  such  as  Fishman  (1986)  or  Karger  (1996).  Other  algorithms 
determine  some  or  all  of  the  coefficients  of  the  generator,  such  as  Colbourn,  Debroni,  and  Myrvold  (1988). 
These  problems  arc  related,  but  are  not  equivalent  especially  in  terms  of  measuring  error.  For  example,  a 
fully-polynomial  randomized  approximation  scheme  (fpras)  for  evaluating  r(xo)  where  xo  is  held  constant 
does  not  yield  a  fpras  for  the  coefficients  of  p.  The  algorithm  of  Beichl,  Cloteaux,  and  Sullivan  (2010), 
which  we  denote  TOP-DOWN,  is  a  sequential  importance  sampling  (SIS)  algorithm.  See  Liu  (2001)  for 
an  explanation  of  SIS. 

TOP-DOWN  ALGORITHM 

Set  mf-  1 . 

For  k  =  0,...,K  do 

Set  Nk  4—  m/k\ 

Determine  D ,  the  set  of  available  non-bridges  of  G . 

Choose  an  edge  e  £  D  uniformly . 

Set  m  4—  m\D\ . 

Set  G<-G-e 

Endf or 

This  pseudocode  is  for  one  estimate.  Multiple  samples  can  be  averaged  to  yield  an  estimate  Nk  where 
we  use  the  bar  to  indicate  a  mean.  We  draw  as  many  samples  as  necessary  until  a  reasonable  accuracy  is 
achieved.  We  may  describe  this  as  a  top-down  algorithm,  as  it  starts  with  the  original  graph  and  removes 
edges  until  reaching  a  tree.  This  algorithm  as  presented  has  running  time  0(\E\2).  However,  as  described  in 
Harris,  Sullivan,  and  Beichl  (2011),  the  running  time  can  be  reduced  to  approximately  0(E\og{V)a(V)), 
where  a()  is  the  inverse  Ackermann  function.  To  simplify  the  notation  in  this  paper,  we  will  abuse  notation 
so  that  E  and  V  can  mean  \E\  and  |V|.  We  also  let  t(G)  represent  the  number  of  spanning  trees  of  a  graph 
G. 

While  this  top-down  algorithm  is  very  effective  at  estimating  the  low-order  coefficients  of  the  reliability 
generator,  its  relative  variance  tends  to  increase  on  the  high-order  coefficients.  In  many  physical  situations, 
only  the  low-order  coefficients  of  the  reliability  generator  arc  important.  For  example,  we  would  likely  be 
interested  in  knowing  when  the  network  has  a  non-negligible  chance  of  staying  connected  after  a  series 
of  disruptions.  Later  coefficients,  corresponding  to  a  large  number  of  broken  links  and  an  exponentially 
small  probability  of  connectivity,  would  probably  not  be  physically  relevant.  However,  in  cases  where  the 
high-order  coefficients  arc  desired,  Colbourn,  Debroni,  and  Myrvold  (1988)  proposed  an  alternative  SIS 
algorithm.  This  algorithm,  which  we  denote  BOTTOM-UP-UNIFORM,  can  be  viewed  as  a  “bottom-up” 
algorithm:  it  starts  with  a  spanning  tree  of  the  graph  and  inserts  edges.  We  sketch  a  simplified  version  of 
the  algorithm  described  in  Colbourn,  Debroni,  and  Myrvold  (1988).  Their  algorithm  has  some  additional 
optimizations  which  will  not  be  relevant  in  this  paper.) 

BOTTOM-UP-UNIFORM  ALGORITHM 

Choose  a  spanning  tree  H  of  G  uniformly  at  random. 

For  k  =  K, . . .  ,0  do 

Set4^t(G)|py 

Choose  an  edge,  e,  of  G  —  H  uniformly  at  random. 

Set  H  4—  H Ue 

Endf or 

This  algorithm  tends  to  have  low  relative  variance  on  high  order  coefficients.  Choosing  the  spanning 
tree  uniformly  at  random  can  be  done  with  Wilson’s  algorithm  (Wilson  1996). 
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Some  reliability  algorithms  return  a  single  estimate  p(x)  of  the  reliability  generator.  In  this  paper  we 
arc  estimating  the  entire  polynomial,  not  simply  evaluating  the  polynomial  at  particular  value  of  x,  and  so 
our  estimate  actually  consists  of  all  K  coefficients. 

In  this  paper,  we  will  analyze  and  improve  the  BOTTOM-UP-UNIFORM.  We  call  the  new  method 
BOTTOM-UP-NEW.  We  make  two  key  improvements.  First,  the  BOTTOM-UP-UNIFORM  algorithm 
requires  computing  the  number  of  spanning  trees  of  certain  subgraphs  of  the  graph  G.  Computing  this 
directly  is  very  expensive,  costing  0(EV 3)  runtime.  We  will  describe  a  more  efficient  method,  which 
speeds  up  BOTTOM-UP-UNIFORM  leading  to  a  runtime  of  G(£2polylog(£.  £)).  It  should  be  noted  that 
Colbourn,  Myrvold,  and  Neufeld  (1996)  uses  a  different  but  similar  linear  algebra  technique  for  unranking 
arboresences. 

The  second  improvement  is  in  the  choice  of  edge  e.  Algorithm  BOTTOM-UP-UNIFORM  selects 
among  the  possible  edges  uniformly  (all  edges  arc  equally  likely).  We  will  describe  an  importance 
sampling  technique  for  selecting  the  edges,  which  greatly  improves  the  accuracy. 

In  a  later  paper  we  will  analyze  the  variance  of  BOTTOM-UP-NEW,  from  a  theoretical  point  of  view 
(including  worst-case  variance).  Here  we  concentrate  on  practical  problems.  In  a  later  paper  we  will  also 
describe  how  to  combine  the  TOP-DOWN  with  BOTTOM-UP-NEW,  in  effect  achieving  the  best  of  both 
worlds  —  low  variance  on  both  high-order  and  low-order  coefficients. 

2  ESTIMATING  HIGH-ORDER  COEFFICIENTS 

As  shown  in  Beichl,  Cloteaux,  and  Sullivan  (2010),  Algorithm  TOP-DOWN  has  very  good  accuracy  on  the 
low-order  coefficients  of  the  reliability  generator.  However  the  algorithm  has  increasing  relative  variance 
for  the  high-order  coefficients.  This  is  a  consequence  of  its  top-down  design.  This  is  unfortunate  because 
the  highest-order  coefficient  of  the  reliability  generator,  which  corresponds  to  the  number  of  subtrees  of 
G,  can  be  computed  exactly  in  polynomial  time  using  the  well-known  Kirchoff  formula.  See,  for  example 
Harris,  Hirst,  and  Mossinghoff  (2008).  Let  Ac  be  the  adjacency  matrix  of  G,  and  let  D  be  a  diagonal 
matrix  whose  z'th  entry  is  the  degree  of  vertex  v;.  Recall  that  t(G)  is  the  number  of  spanning  trees  of  G. 
We  have  the  identity 

t(G)  =  detL, 

where  L  =  (D—Aq)u  denotes  the  minor  of  D—Aq  obtained  by  removing  the  first  row  and  column.  This 
can  be  evaluated  in  0(V3)  work. 

Next,  consider  how  one  might  seek  to  estimate  the  penultimate  coefficient  of  the  reliability  generator. 
This  coefficient  counts  the  number  of  graphs  with  V  (=  (V  —  1)  + 1)  edges.  To  obtain  such  a  graph,  one 
might  take  a  spanning  tree  T  and  add  any  of  the  edges  of  other  K  edges  of  G.  However,  this  would  count 
graphs  multiple  times;  each  such  graph  T  Ue  would  be  counted  with  multiplicity  z(T  Lie).  Accounting  for 
this  factor,  one  has  that 

penultimate  coefficient  of  pc  =  z{G){E  —  V  +  l)E[l/T(rUe)] 

where  the  expectation,  E,  is  taken  over  all  spanning  trees  T  and  all  edges  e  £  G  —  T . 

To  estimate  this  quantity  unbiasedly,  one  might  sample  spanning  trees  T  and  edges  e  G  G  —  T ,  compute 
the  resulting  z(T  Ue),  and  extract  a  sample  mean  of  I  jx(T  LJ  e). 

In  general,  coefficient  K  —  k  of  the  reliability  generator  is  given  by  the  formula 

T(G)E[  £  l/T(7’Ue1-'Uejt)] 

e\,...,ek 

where  T  is  chosen  uniformly  among  spanning  trees  of  G,  e\,.  .  .  ,  e/(  arc  distinct  edges  of  G  —  T. 

This  formula  suggests  that  there  arc  two  search  spaces  to  cover:  the  space  of  spanning  trees  T  and  the 
space  of  edges  e\, . . .  ,e*.  Our  algorithm  will  handle  the  first  using  simple  Monte  Carlo  sampling,  and  will 
handle  the  second  using  Sequential  Importance  Sampling.  This  will  be  a  bottom-up  algorithm:  we  start 
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from  a  spanning  tree  and  add  edges  to  it  until  recovering  the  original  graph.  We  will  denote  this  Algorithm 
BOTTOM-UP-UNIFORM: 

Algorithm  BOTTOM-UP-UNIFORM 


Precompute  z(G) 

Choose  a  spanning  tree  H  of  G  uniformly  at  random 
algorithm) 

For  k  =  K, . . .  ,0  do 


Set  Nk 


r(G)( 


z(H) 

Choose  an  edge  of  G 
Set  H  4—  H Ue 
Endf or 


H  uniformly  at  random. 


(e.g. 


with  Wilson' s 


Algorithm  BOTTOM-UP-UNIFORM  is  equivalent  to  the  algorithm  described  in  Colbourn,  Debroni, 
and  Myrvold  (1988).  In  this  algorithm  as  presented,  as  in  many  SIS  algorithms,  variance  is  introduced 
when  different  choices  for  e  at  each  step  lead  to  estimates  of  divergent  size.  (Variance  is  also  introduced 
with  choice  of  spanning  tree,  but  we  will  not  use  the  importance  to  deal  with  this.)  In  general,  in  SIS 
algorithms,  it  is  best  to  choose  an  importance  function  which  equalizes,  as  far  as  possible,  the  contributions 
from  each  choice.  In  our  case,  this  is  not  really  possible,  for  two  reasons.  First,  we  cannot  compute  the 
exact  contribution  —  that  is  precisely  the  quantity  we  arc  trying  to  estimate  in  the  first  place.  Second,  we 
are  seeking  to  estimate  the  entire  reliability  generator,  and  so  our  selection  contributes  to  multiple  estimated 
coefficients,  which  we  cannot  equalize  simultaneously. 

In  our  case,  we  will  use  a  simple  importance  function  so  that  the  contribution  from  each  edge  e  with 
respect  to  the  next  estimate  Nk  is  equalized.  Letting  Pie)  denote  the  probability  of  inserting  edge  e,  the 
choice  of  P(e )  that  achieves  this  is 


P(e)oc 


1 

T(ffUe) 


which  yields  the  following  algorithm: 

Algorithm  BOTTOM-UP-NEW 


Precompute  Z (G) . 

Set  mi—  1,  the  running  total  of  all  importance  functions  seen  so  far. 
Choose  a  spanning  tree  H  of  G  uniformly  at  random  (e.g.  withWilson's 
algorithm) . 

For  k  =  K, . . .  ,0  do 
Set  Nk^mz(G)(Kk) 

Set 

S=  £  z(H)/z(HUd) 

deG~H 

Update 

m  i—  mS  /  k 

Choose  an  edge  e^G  —  H  with  probability  P{e)  proportional  to  l/z(HUe) . 
Set  H  i—  HUe 
Endf  or 


In  practice.  Algorithm  BOTTOM-UP-NEW  has  much  better  variance  than  Algorithm  BOTTOM-UP- 
UNIFORM  on  the  high-order  coefficients. 
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The  results  of  Algorithm  BOTTOM-UP-NEW  arc  almost  the  opposite  of  Algorithm  TOP-DOWN: 
high-order  coefficients  arc  computed  very  accurately  (in  fact,  the  highest-order  coefficient  is  computed 
exactly),  while  the  relative  variance  increases  for  the  low-order  coefficients. 

3  COMPUTING  X  EFFICIENTLY 

Algorithms  BOTTOM-UP-UNIFORM  and  BOTTOM-UP-NEW  appeal-  to  be  very  expensive  because  they 
requires  extensive  computation  of  the  x  function.  Directly  computing  x  using  the  Kirchhoff  determinant 
formula  costs  0(V3)  work  and  0(V2)  memory  (for  Gaussian  elimination).  It  would  appear  that  at  each 
iteration,  we  must  compute  x(GUe)  for  each  edge;  and  altogether  there  are  K  iterations.  In  all,  this  would 
appeal-  to  cost  0(E2V 3)  work  and  0(V2)  memory  for  Algorithm  BOTTOM-UP-NEW. 

Instead  of  recomputing  x(G Lie)  at  each  stage,  we  will  instead  keep  track  of  the  value 

1(e)  =  x(HUe)/x(H) 

at  each  stage  of  the  algorithm,  for  each  edge  e  G  G  —  H. 

Initially,  when  H  is  a  spanning  tree,  we  set 

X(i,j)  tree-distance  (ij)  +  1 


for  each  edge  (i,j). 

Next,  when  we  update  H  by  inserting  an  edge  e,  we  must  update  A.  Let  d  be  another  edge  in  H .  We 
need  to  compute  the  updated  X'(d )  =  x(H  U eUd) /x(H L) e). 

We  will  find  the  following  notation  convenient.  If  e  =  (i,j)  is  any  edge  in  G,  we  define  the  column 
vector  8e  to  be  the  vector  which  is  +1  in  coordinate  i,  is  —1  in  coordinate  j,  and  is  zero  elsewhere.  Observe 
that  when  edge  e  is  added  to  H ,  the  matrix  L  changes  by  8e8j : 

1} !  e  —  l j  H  T  8e8e  . 

Now  let  us  examine  how  to  recompute  A': 


A  '(d) 


x(HUeUd)/x(HUe) 

det  (Lh  +  8d8j  +  8e8j )/ det  (LH  +  8e8j ) 


T'~x8d8T 


det  [I  +  (Lh  +  8e  8e  j  -  od  od 

1  +  8d  (Lh  +  8eSj )  1 8d 

j~l 8  8T 

i  +  STALH'-L«  b±L«)S„ 


I  -  8j  Ln'  8„ 


1  +  5JW- 


lilt 


1  -8T\ 


-)8d 


where  u  =  LHl8e 


1  +8dLHl8d  - 


(gJ“)2 

1  —  8J  u 


A  (d) 


(gj»)2 

m  ' 


This  leads  to  the  following  technique  for  updating  A.  Whenever  we  add  edge  e  to  H ,  we  use  the  conjugate 
gradient  method,  a  sparse  matrix  technique,  to  compute  u  =  Ln]  8,,.  We  then  update  each  A  (d) 

X(d)^X(d)-(8juf/X(e) 
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In  total,  this  technique  allows  us  to  reduce  the  memory  requirement  to  0(E)  and  the  complexity  to 
0(EV).  Furthermore,  in  practice  the  conjugate  gradient  method  needs  far  less  than  0(V)  iterations  to 
compute  an  adequate  approximation,  so  the  work  requirement  per  iteration  is  more  like  <9(f?polylog(£',  e)) 
where  e  is  the  desired  accuracy  after  averaging  multiple  samples  together. 

To  make  this  completely  transparent,  we  present  Algorithm  BOTTOM-UP-NEW  explicitly  describing 
how  to  use  and  update  X: 

Algorithm  BOTTOM-UP-NEW  (with  linear  algebra) 

Precompute  T (G) . 

Set  mt-  1,  the  running  total  of  all  importance  functions  seen  so  far. 
Choose  a  spanning  tree  H  of  G  uniformly  at  random  (e.g.  with  Wilson' s 
algorithm) . 

For  each  edge  (i,j)  G  G  —  H,  set 

X(i,j)  <—  tree-distance  (i,j)  +  1 

For  k  =  K, . . .  ,0  do 

Set  pk  A-  mz(G)(f) 

Set  S  =  Yjd€G~H  1  /A(d)  and  update  m  mS/k 

Choose  an  edge  e  £G  —  H  with  probability  P(e)  =  (1  IX(e))/S . 

Using  conjugate  gradient  method,  compute  u~Ljf8e 

For  each  edge  d  €  G  —  H  update 

X(d)^X(d)-(8ju)2/ X(e) 

Endf or 

Set  H  <—  HUe 
Endf or 

This  is  mathematically  equivalent  to  the  previous  Algorithm  BOTTOM-UP-NEW.  In  total,  these  tech¬ 
niques  reduce  the  work  requirement  of  both  Algorithm  BOTTOM-UP-NEW  and  Algorithm  BOTTOM- 
UP-UNIFORM  to  a  maximum  of  0(E2V )  and  memory  requirement  to  0(E)',  and  in  practice  much  less 
than  this  if  we  accept  a  numerical  approximation  (more  like  <9(£'2polylog(£’,£))). 

Also,  for  large  values  of  k,  the  matrix  L  is  extremely  sparse.  In  particular,  it  has  many  density-two 
rows  (corresponding  to  degree-one  vertices  in  H ).  For  inverting  a  matrix  of  this  form,  a  simple  application 
of  conjugate  gradient  is  not  best.  A  better  strategy  is  to  use  a  partial  Gaussian  elimination,  eliminating  all 
degree-two  rows,  and  performing  a  preconditioned  conjugate  gradient  on  the  resulting  smaller  matrix.  We 
have  found  that  the  inverse  diagonal  preconditioner  works  well. 

Another  trick  to  speed  up  this  step  is  to  note  that  the  Kirchoff  formula  works  for  any  minor  of  the 
matrix  D  —  A.  To  reduce  the  density  of  the  resulting  matrix  L,  we  chose  to  omit  the  row  and  column  i 
corresponding  to  the  highest-degree  vertex  (after  pruning  away  degree-one  vertices). 

In  practice,  these  tricks  can  speed  up  Algorithm  BOTTOM-UP-NEW  by  a  factor  of  four  or  more  for 
large  k.  For  small  k,  the  matrix  L  is  not  so  sparse  and  these  techniques  gain  little. 

4  RESULTS 

Let  us  turn  to  some  sample  graphs  to  assess  the  performance  of  these  algorithms.  Our  first  test  case 
is  an  Erdos-Renyi  graph  with  100  vertices  and  150  edges.  Figure  1  depicts  the  relative  variance  of  the 
three  Algorithms:  Algorithm  TOP-DOWN,  Algorithm  BOTTOM-UP-NEW  and  Algorithm  BOTTOM-UP- 
UNIFORM. 

As  we  have  discussed  earlier,  the  relative  variance  of  Algorithm  TOP-DOWN  increases  for  large  k 
and  the  relative  variance  of  Algorithm  BOTTOM-UP-NEW  increases  for  small  k.  Algorithm  BOTTOM- 
UP-UNIFORM,  as  we  have  discussed,  has  better  relative  variance  than  Algorithm  BOTTOM-UP-NEW  on 
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Figure  1:  Relative  variance  of  various  algorithms  on  Erdos  Renyi  graph,  100  vertices  250  edges.  Key: 
Solid  =  TOP-DOWN  algorithm.  Dashed  =  BOTTOM-UP-UNIFORM  (uniform  distribution),  Dotted  = 
BOTTOM-UP-NEW  (new  importance  sampling). 

low-order  coefficients  and  worse  variance  on  high-order  coefficients.  Since  Algorithm  TOP-DOWN  also 
has  low  relative  variance  on  these  low-order  coefficients.  Algorithm  BOTTOM-UP-UNIFORM  seems  to 
be  always  dominated  by  either  Algorithm  TOP-DOWN  or  Algorithm  BOTTOM-UP-NEW. 

Test  Case  2  comes  from  a  Delaunay  triangulation  on  100  points,  yielding  a  graph  with  100  vertices  and 
286  edges.  Figure  2  depicts  the  relative  variance  of  Algorithms  TOP-DOWN,  BOTTOM-UP-UNIFORM, 
and  BOTTOM-UP-NEW. 

For  any  coefficient  k,  either  TOP-DOWN  or  BOTTOM-UP-NEW  has  lower  variance  than  BOTTOM- 
UP-UNIFORM.  Furthermore,  either  TOP-DOWN  or  BOTTOM-UP-NEW  is  at  least  as  fast  as  BOTTOM- 
UP-UNIFORM. 

5  CONCLUSION 

We  have  described  a  new  bottom-up  algorithm  for  estimating  the  reliability  polynomial  of  a  connected 
graph.  This  new  algorithm  achieves  substantially  better  speed  and  accuracy  compared  to  the  algorithm 
of  Colbourn,  Debroni,  and  Myrvold  (1988).  This  improved  accuracy  is  especially  beneficial  because 
the  top-down  algorithm  of  Beichl,  Cloteaux,  and  Sullivan  (2010)  has  poor  accuracy  in  the  high-order 
coefficients. 
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Figure  2:  Relative  Variance  of  various  algorithms  on  Delanuay  data.  Key:  Solid  =  TOP-DOWN  algo¬ 
rithm,  Dashed  =  BOTTOM-UP-UNIFORM  (uniform  distribution).  Dotted  =  BOTTOM-UP-NEW  (new 

distribution). 
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